#include "imgpro.h"
using namespace cv;
using namespace std;
void MedianFilter(cv::Mat M,cv::Mat M1, int window)
{
    int row = M.rows;
    int col = M.cols;
    int step0 = M.step[0];
    int step1 = M.step[1];
    for (int i = 0; i < row; i++)
        for (int j = 0; j < col; j++)
        {
            Vec3b *a = new Vec3b[window*window];
            int k = 0;
            for (int ii = -window / 2; ii <= window / 2; ++ii)
                for (int jj = -window / 2; jj <= window / 2; ++jj){
                    if ((i + ii < 0) || (j + jj < 0) || (i + ii >= row) || (j + jj >= col))
                        continue;
                    a[k][0] = (M.data + (i + ii) * step0 + (j + jj) * step1)[0];
                    a[k][1] = (M.data + (i + ii) * step0 + (j + jj) * step1)[1];
                    a[k++][2] = (M.data + (i + ii) * step0 + (j + jj) * step1)[2];

                }
            for (int m = 0; m < 3; ++m)
            {
                int count[256] = { 0 };
                for (int i = 0; i < k; ++i)
                    count[a[i][m]]++;
                int medium = 1;
                for (; medium < 256; ++medium)
                {
                    if (count[medium - 1] >= k / 2 + 1)
                        break;
                    count[medium] += count[medium - 1];
                }
                for (int jj = 0; jj < k; ++jj)
                {
                    if (a[jj][m] + 1 == medium)
                    {
                        (M1.data + i  * step0 + j  * step1)[m] = a[jj][m];
                        break;
                    }
                }
            }
        }
}
void GaussianMat(double*a, int n, int delta)
{
    double sum = 0; int k = 0;
    for (int i = -n / 2; i <= n / 2; ++i)
        for (int j = -n / 2; j <= n / 2; ++j)
        {
            a[k] = pow(M_E, -0.5 * (i*i + j*j) / (delta*delta))*0.5 / (M_PI*delta*delta);
            sum += a[k++];
        }
    for (int i = 0; i < k; ++i)
        a[i] = a[i] / sum;

}



void GaussianBlur(Mat M, Mat M1, int window,int delta)
{
    int row = M.rows;
    int col = M.cols;
    int step0 = M.step[0];
    int step1 = M.step[1];
    double *gaussianMat = new double[window*window];
    GaussianMat(gaussianMat, window, delta);
    for (int i = 0; i < row; i++)
        for (int j = 0; j < col; j++)
        {
            double a[3] = { 0 };
            int k = 0;
            for (int ii = -window / 2; ii <= window / 2; ++ii)
                for (int jj = -window / 2; jj <= window / 2; ++jj){
                    if ((i + ii < 0) || (j + jj < 0) || (i + ii >= row) || (j + jj >= col))
                        continue;
                    a[0] += (M.data + (i + ii) * step0 + (j + jj) * step1)[0] * gaussianMat[k];
                    a[1] += (M.data + (i + ii) * step0 + (j + jj) * step1)[1] * gaussianMat[k];
                    a[2] += (M.data + (i + ii) * step0 + (j + jj) * step1)[2] * gaussianMat[k++];
                }
            if (k != 0)
            {
                (M1.data + i  * step0 + j  * step1)[0] = a[0];
                (M1.data + i  * step0 + j  * step1)[1] = a[1];
                (M1.data + i  * step0 + j  * step1)[2] = a[2];
            }
        }
}
void Filter(Mat M,Mat M1, int window)
{
    int row = M.rows;
    int col = M.cols;
    int step0 = M.step[0];
    int step1 = M.step[1];
    for (int i = 0; i < row; i++)
        for (int j = 0; j < col; j++)
        {
            int a[3] = { 0 };
            int k = 0;
            for (int ii = -window / 2; ii <= window / 2; ++ii)
                for (int jj = -window / 2; jj <= window / 2; ++jj){
                    if ((i + ii < 0) || (j + jj < 0) || (i + ii >= row) || (j + jj >= col))
                        continue;
                    a[0] += (M.data + (i + ii) * step0 + (j + jj) * step1)[0];
                    a[1] += (M.data + (i + ii) * step0 + (j + jj) * step1)[1];
                    a[2] += (M.data + (i + ii) * step0 + (j + jj) * step1)[2];
                    ++k;
                }
            if (k != 0)
            {
                (M1.data + i  * step0 + j  * step1)[0] = a[0] / k;
                (M1.data + i  * step0 + j  * step1)[1] = a[1] / k;
                (M1.data + i  * step0 + j  * step1)[2] = a[2] / k;
            }
        }
}
void SaturateProc(Mat M,Mat M1, int window)
{
    int row = M.rows;
    int col = M.cols;
    int step0 = M.step[0];
    int step1 = M.step[1];
    cv::cvtColor(M, M1, CV_BGR2HSV);
    for (int i = 0; i < row; i++)
        for (int j = 0; j < col; j++)
        {
            if((M1.data + i * step0 + j * step1)[1]+window<0)
                (M1.data + i * step0 + j * step1)[1]=0;
            else if((M1.data + i * step0 + j * step1)[1]+window>255)
                (M1.data + i * step0 + j * step1)[1]=255;
            else
                (M1.data + i * step0 + j * step1)[1]+=window;

        }
    cv::cvtColor(M1, M1, CV_HSV2BGR);

}

void ghostFilter(Mat M,Mat M1, int window)
{
    int row = 512;
    int col = 512;
    for (int i = 0; i < col; i++)
        for (int j = 0; j < row; j++)
        {
            Vec3b *a = new Vec3b[window*window];
            int k = 0;
            for (int ii = -window / 2; ii <= window / 2; ++ii)
                for (int jj = -window / 2; jj <= window / 2; ++jj){
                    if ((i + ii < 0) || (j + jj < 0) || (i + ii >= col) || (j + jj >= row))
                        continue;
                    a[k][0] = M.at<Vec3b>(i + ii, j + jj)[0];
                    a[k][1] = M.at<Vec3b>(i + ii, j + jj)[1];
                    a[k++][2] = M.at<Vec3b>(i + ii, j + jj)[2];
                }
            for (int m = 0; m < 3; ++m)
            {
                int count[256] = { 0 };
                for (int i = 0; i < k; ++i)
                    count[a[i][m]]++;
                int medium = 1;
                for (; medium < 256; ++medium)
                {
                    if (count[medium - 1] >= k / 2)
                        break;
                    count[medium] += count[medium - 1];
                }
                for (int jj = 0; jj < k; ++jj)
                {
                    if (a[jj][m] + 1 == medium)
                    {
                        M1.at<Vec3b>(i, j)[m] = a[jj][m];
                        break;
                    }
                }
            }
        }
}
void sharpen(Mat img, Mat result)
{
    result.create(img.size(), img.type());
    for (int row = 1; row < img.rows - 1; row++)
    {
        const uchar* previous = img.ptr<const uchar>(row - 1);
        const uchar* current = img.ptr<const uchar>(row);
        const uchar* next = img.ptr<const uchar>(row + 1);
        uchar *output = result.ptr<uchar>(row);
        int ch = img.channels();
        int starts = ch;
        int ends = (img.cols - 1) * ch;
        for (int col = starts; col < ends; col++)
        {
            *output++ = saturate_cast<uchar>(9	 * current[col] - current[col - ch]
                    - current[col + ch] - previous[col] - next[col]-previous[col-ch]-previous[col+ch]
                    - next[col-ch]-next[col+ch]);
        }
    }
}


void  mapping(warp *w, int x, int y, double *u, double *v)
{
    double fu, fv;
    fu = x;
    fv = y;
    double dx = fu - w->origCenterX;
    double dy = fv - w->origCenterY;
    double rsq = dx*dx + dy*dy;
    if (dx>0.0 - w->radius
        && dx<w->radius
        && dy>0.0 - w->radius
        && dy<w->radius
        && rsq<  w->radius*w->radius
        )
    {
        double cmx = w->nowCenterX - w->origCenterX;
        double cmy = w->nowCenterY - w->origCenterY;
        double msq = cmx*cmx + cmy*cmy;
        double edge_dist = w->radius*w->radius - rsq;
        double a = edge_dist / (edge_dist + msq);
        a *= a;
        fu -= a*cmx;
        fv -= a*cmy;
    }
    *u = fu;
    *v = fv;
}
void imageWarp(warp * w)
{
    IplImage temp = (IplImage)w->img;
    IplImage *img = &temp;
    temp = (IplImage)w->dst;
    IplImage *dst = &temp;
    for (int i = 0; i<img->height; i++)
    {
        uchar *ptr = (uchar*)(img->imageData + i*img->widthStep);
        uchar *ptrWarp = (uchar *)(dst->imageData + i*dst->widthStep);
        for (int j = 0; j<img->width; j++)
        {
            double dx = j - w->origCenterX;
            double dy = i - w->origCenterY;
            if (dx>(0.0 - w->radius)
                && dx<(w->radius)
                && dy>(0.0 - w->radius)
                && dy<w->radius
                && (dx*dx + dy*dy) < w->radius *w->radius
                )
            {

                double u, v;
                mapping(w, j, i, &u, &v);
                if (u>img->width)
                    u = img->width;
                if (v>img->height)
                    v = img->height;
                if (u<0)
                    u = 0;
                if (v<0)
                    v = 0;
                uchar *ptr1 = (uchar *)(IplImage(w->img).imageData + (int)v*img->widthStep);
                int iu = (int)u;
                ptrWarp[3 * j + 0] = ptr1[3 * iu + 0];
                ptrWarp[3 * j + 1] = ptr1[3 * iu + 1];
                ptrWarp[3 * j + 2] = ptr1[3 * iu + 2];

            }
            else
            {
                ptrWarp[3 * j + 0] = ptr[3 * j + 0];
                ptrWarp[3 * j + 1] = ptr[3 * j + 1];
                ptrWarp[3 * j + 2] = ptr[3 * j + 2];
            }
        }
    }
    //cvCopy(dst, &IplImage(w->dst));
    temp = (IplImage)w->img;
    cvCopy(dst, &temp);
}

bool R1(int R, int G, int B) {
    bool e1 = (R>95) && (G>40) && (B>20) && ((max(R, max(G, B)) - min(R, min(G, B)))>15) && (abs(R - G)>15) && (R>G) && (R>B);
    bool e2 = (R>220) && (G>210) && (B>170) && (abs(R - G) <= 15) && (R>B) && (G>B);
    return (e1 || e2);
}

bool R2(float Y, float Cr, float Cb) {
    bool e3 = Cr <= 1.5862*Cb + 20;
    bool e4 = Cr >= 0.3448*Cb + 76.2069;
    bool e5 = Cr >= -4.5652*Cb + 234.5652;
    bool e6 = Cr <= -1.15*Cb + 301.75;
    bool e7 = Cr <= -2.2857*Cb + 432.85;
    return e3 && e4 && e5 && e6 && e7;
}

bool R3(float H, float S, float V) {
    return (H<25) || (H > 230);
}

Mat GetSkin(Mat const &src) {
    // allocate the result matrix
    Mat dst;
    src.copyTo(dst);

    Vec3b cwhite = Vec3b::all(255);
    Vec3b cblack = Vec3b::all(0);

    Mat src_ycrcb, src_hsv;
    // OpenCV scales the YCrCb components, so that they
    // cover the whole value range of [0,255], so there's
    // no need to scale the values:
    cvtColor(src, src_ycrcb, CV_BGR2YCrCb);
    // OpenCV scales the Hue Channel to [0,180] for
    // 8bit images, so make sure we are operating on
    // the full spectrum from [0,360] by using floating
    // point precision:
    src.convertTo(src_hsv, CV_32FC3);
    cvtColor(src_hsv, src_hsv, CV_BGR2HSV);
    // Now scale the values between [0,255]:
    normalize(src_hsv, src_hsv, 0.0, 255.0, NORM_MINMAX, CV_32FC3);

    for (int i = 0; i < src.rows; i++) {
        for (int j = 0; j < src.cols; j++) {

            Vec3b pix_bgr = src.ptr<Vec3b>(i)[j];
            int B = pix_bgr.val[0];
            int G = pix_bgr.val[1];
            int R = pix_bgr.val[2];
            // apply rgb rule
            bool a = R1(R, G, B);

            Vec3b pix_ycrcb = src_ycrcb.ptr<Vec3b>(i)[j];
            int Y = pix_ycrcb.val[0];
            int Cr = pix_ycrcb.val[1];
            int Cb = pix_ycrcb.val[2];
            // apply ycrcb rule
            bool b = R2(Y, Cr, Cb);

            Vec3f pix_hsv = src_hsv.ptr<Vec3f>(i)[j];
            float H = pix_hsv.val[0];
            float S = pix_hsv.val[1];
            float V = pix_hsv.val[2];
            // apply hsv rule
            bool c = R3(H, S, V);

            if (!(a&&b&&c))
                dst.ptr<Vec3b>(i)[j] = cblack;
        }
    }
    return dst;
}


void _BilateralFilter(Mat M, Mat M1,int deltad ,int deltar, int window)
{
    int row = M.rows;
    int col = M.cols;
    int step0 = M.step[0];
    int step1 = M.step[1];
    int pixel0[3]={0};
    for (int i = 0; i < row; i++)
        for (int j = 0; j < col; j++)
        {
            double sum[3] = { 0 } ,b[3] = {0},bsum[3]={0} , gaussian=0;
            int pixeln[3]={0};
            pixel0[0] = (M.data + i * step0 + j * step1)[0];
            pixel0[1] = (M.data + i * step0 + j * step1)[1];
            pixel0[2] = (M.data + i * step0 + j * step1)[2];
            for (int ii = -window / 2; ii <= window / 2; ++ii)
            {
                if ((i + ii < 0) || (i + ii >= row))
                    continue;
                gaussian = pow(M_E, -0.5*(ii*ii)/(deltad*deltad));
                for(int k=0;k<3;++k)
                {
                    pixeln[k] = (M.data + (i + ii) * step0 + (j) * step1)[k];
                    b[k] = pow(M_E, -0.5*(pixeln[k]-pixel0[k])*(pixeln[k]-pixel0[k])/(deltar*deltar))*gaussian;
                    bsum[k] += b[k];
                    sum[k] += pixeln[k]*b[k];
                }
            }
            for(int k=0;k<3;++k)
                (M1.data + i  * step0 + j  * step1)[k] = sum[k]/bsum[k];


            sum[3] = { 0 } ,b[3] = {0},bsum[3]={0} , gaussian=0;
            for (int jj = -window / 2; jj <= window / 2; ++jj)
            {
                if ((j + jj < 0) || (j + jj >= col))
                    continue;
                gaussian = pow(M_E, -0.5*(jj*jj)/(deltad*deltad));
                for(int k=0;k<3;++k)
                {
                    pixeln[k] = (M.data + i * step0 + (j+jj) * step1)[k];
                    b[k] = pow(M_E, -0.5*(pixeln[k]-pixel0[k])*(pixeln[k]-pixel0[k])/(deltar*deltar))*gaussian;
                    bsum[k] += b[k];
                    sum[k] += pixeln[k]*b[k];
                }
            }
            for(int k=0;k<3;++k)
                (M1.data + i  * step0 + j  * step1)[k] = sum[k]/bsum[k];
        }
}

void BilateralFilter(Mat M, Mat M1,int deltad ,int deltar, int window)
{
    Mat s=GetSkin(M);
    Mat skin;
    skin.create(M.size(), M.type());
    M.copyTo(M1);
    _BilateralFilter( s, skin,deltad ,deltar, window);
    int rows=M1.rows;
    int cols=M1.cols;
    int step0 = M1.step[0];
    int step1 = M1.step[1];
    int a0 = 0, a1 = 0, a2 = 0;
    for (int i = 0; i < rows; ++i)
        for (int j = 0; j < cols; ++j)
        {
               a0 = (skin.data + i * step0 + j * step1)[0];
               a1 = (skin.data + i * step0 + j * step1)[1];
               a2 = (skin.data + i * step0 + j * step1)[2];
               if (a0 + a1 + a2 < 150)
                   continue;
               else
               {
                   (M1.data + i * step0 + j * step1)[0] = a0;
                   (M1.data + i * step0 + j * step1)[1] = a1;
                   (M1.data + i * step0 + j * step1)[2] = a2;

               }
         }
}


void WhiteShit(Mat M, Mat M1, int alpha)
{
    M.copyTo(M1);
    alpha = -alpha;
    Mat grey;
    cvtColor(M, grey, CV_BGR2GRAY);
    int x, y;
    uchar*Src, *Dst, *gg;
    for (y = 0; y < M1.rows; ++y)
    {
        Src = M.data + y*M.step[0];
        Dst = M1.data + y*M1.step[0];
        gg = grey.data + y*grey.step[0];
        for (x = 0; x < M.cols; ++x)
        {
            if ((int)(gg[0] + alpha) >= 0&&(int)(gg[0] + alpha) <= 255)
                gg[0] += alpha;
            else if((int)(gg[0] + alpha) <= 0)
                gg[0] = 0;
            else
                gg[0] = 255;

            if (Src[0]>20 && Src[1] > 40 && Src[2] > 50)
            {
                Dst[0] = (Src[0] * (255 - gg[0]) + 255 * gg[0]) / 255;
                Dst[1] = (Src[1] * (255 - gg[0]) + 255 * gg[0]) / 255;
                Dst[2] = (Src[2] * (255 - gg[0]) + 255 * gg[0]) / 255;
            }
            Src += 3;
            Dst += 3;
            gg++;
        }
    }
}

void LomoFilter(Mat M,Mat M1,int window)
{
    float mSize = 0.5;
    int width = M.cols;
    int heigh = M.rows;
    int centerX = width >> 1;
    int centerY = heigh >> 1;

    int maxV = centerX*centerX + centerY*centerY;
    int minV = (int)(maxV*(1 - mSize));
    int diff = (maxV - minV);
    float ratio = width >heigh ? (float)heigh / (float)width : (float)width / (float)heigh;
    for (int y = 0; y<heigh; y++)
    {
        uchar* imgP = M.ptr<uchar>(y);
        uchar* dstP = M1.ptr<uchar>(y);
        for (int x = 0; x<width; x++)
        {
            int b = imgP[3 * x]*0.85 ;
            int g = imgP[3 * x + 1];
            int r = imgP[3 * x + 2];
            float dx = centerX - x;
            float dy = centerY - y;
            if (width > heigh)
                dx = (dx*ratio);
            else
                dy = (dy*ratio);

            int dstSq = (dx*dx + dy*dy);
            float v = ((float)dstSq / diff) * window;
            r = (int)(r - v);
            g = (int)(g - v);
            b = (int)(b - v);
            r = (r>255 ? 255 : (r<0 ? 0 : r));
            g = (g>255 ? 255 : (g<0 ? 0 : g));
            b = (b>255 ? 255 : (b<0 ? 0 : b));
            dstP[3 * x] = (uchar)b;
            dstP[3 * x + 1] = (uchar)g;
            dstP[3 * x + 2] = (uchar)r;

        }
    }

}

